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Abstract 

The water-graphite interaction potential proposed recently (Gonzalez et al. J. 
Phys. Chem. C 2007, 111, 14862), the three TIPiVP (N = 3, 4, 5) water- water 
interaction models, and basin-hopping global optimization are used to find the likely 
candidates for the global potential energy minima of (H20) n clusters with n < 21 on 
the (OOOl)-surface of graphite and to perform a comparative study of these minima. 
We show that, except for the smaller clusters (n < 6), for which ab-initio results 
are available, the three water-water potential models provide mostly inequivalent 
conformations. While TIP3P seems to favor monolayer water structures for n < 18, 
TIP4P and TIP5P favor bilayer or volume structures for n > 6. These n values 
determine the threshold of dominance of the hydrophobic nature of the water-graphite 
interaction at the nanoscopic scale for these potential models. 



1 Introduction 

The interaction between water and graphite has been the concern of theoretical and ex- 
perimental studies. A deep understanding of the features and properties of this interaction 
is of great interest in technological applications pflElB], environmental sciences [4], and 
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astrophysics [5j, among other fields. The establishment of either the hydrophilic or hy- 
drophobic nature of graphite at nanoscopic scales, which is of particular relevance in those 
applications, must be based on the knowledge of this interaction. 

In a previous publication [6 J (hereafter referred to as I), we have developed a model 
for the water-graphite interaction and found the likely candidates for the global potential 
energy minima of (H 2 0) n clusters with n < 21 on the (OOOl)-surface of graphite. Out 
of this model, we have obtained a rather hydrophobic water-graphite interaction at the 
nanoscopic scale. As a consequence of this property, the water component of the lowest 
graphite-(H20) n minima is quite closely related to low-lying minima of the corresponding 
(H 2 0) n clusters. In about half of the cases the geometrical substructure of the water 
molecules in the graphite- (H20) n global minimum coincides with that of the corresponding 
free water cluster. Exceptions occur when the interaction with graphite induces a change 
in the geometry of the water moiety. Our general conclusions were in agreement with the 
sparse experimental [8] and theoretical data 0, M, HH [HI CE2] . Besides, the structures of 
these minima for 1 < n < 6 coincided with those provided by empirical pEU] and ab initio 
calculations [3l \TT\. 

In our study, the water-water interaction was described by the TIP4P intermolecular 
potential model [1 3J . The related TIP3P potential p[3] was also used to model this in- 
teraction for n < 6. The global minimum structures found for these clusters coincided 
with those of the TIP4P model. However, the observed dependence of the structure of the 
water-graphite global minima on the structure of the corresponding free water clusters and 
the known dependence of the latter on the water-water interaction model for n > 6 antic- 
ipated a dependence of the structure of these larger water-graphite clusters on the form 
chosen to model the water-water interaction. In I, preliminary results with the TIP3P 
model confirmed this prediction. In this article, we will present the concluding results 
from our analysis of this dependence by considering also the TIP5P model [14]. As in I, 
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we will make use of the water-graphite interaction model developed there and the basin- 
hopping method to find the likely candidates for the global potential energy minima of 
graphite-(H20)„ clusters with n < 21 and the TIP3P and TIP5P water-water interaction 
models, and perform a systematic comparison of the cluster structures found with these 
and the TIP4P model. 

This paper is organized as follows. In Section [2] we summarize the relevant details of 
the model developed in I for the water-graphite interaction. In Section [3] we present likely 
candidates for the cluster global potential energy minima together with their association 
and binding energies for both the TIP3P and TIP5P water-water interaction models. We 
shall also compare these global minimum structures with those found in I for the TIP4P 
model. Finally, Section [4] summarizes our conclusions. 

2 Summary of the Potential Energy Function 

The closed-shell electronic structure of both graphite and water makes an empirical ap- 
proach to the potential energy surface (PES) for the water-graphite and water-water in- 
teractions particularly attractive. In I, we wrote the potential energy of a graphite- (H 2 0) n 
cluster as a sum of two contributions 

V = V WW + V wg , (1) 

where V ww is the sum of pairwise water-water interactions, and V wg is the water-graphite 
term. For the water-water interaction, the TIP4P model was the primary choice in I; 
here we will study the performance of the TIP3P and TIP5P potentials. All these models 
describe each water molecule as the same rigid body with two positive charges on the 
hydrogen atoms and either a balancing negative charge at the oxygen atom (TIP3P) or 
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close to the oxygen atom (TIP4P), or two balancing negative charges close to the oxygen 
atom and out of the molecular plane (TIP5P), together with a dispersion-repulsion center 
on the oxygen atom. Hence, V ww is a sum of pairwise additive Coulomb and Lennard- 
Jones terms. We should remind here that the TIPiVP are a family of empirical water- 
water potentials whose parameters have been appropriately set so as to reproduce some 
properties of the liquid water phase at room temperature. Potentials from these family 
have been used in the study of homogeneous water clusters [H2 [EH [TZ1 [[8] , water clusters 
containing metallic cations [El [20], and water-Ceo clusters [21]. 
The water-graphite interaction is written as 

V wg = V dv + V poh (2) 

where V& r is a sum of pairwise dispersion-repulsion terms between the oxygen and the 
carbon atoms. Each of these terms is expressed as a Lennard- Jones potential, whose 
parameters were obtained using the standard Lorentz-Berthelot combination rules from the 
corresponding parameters for the oxygen-oxygen and carbon-carbon interactions in TIPiVP 
water and Steele [22] graphene-graphene potentials, respectively. Specifically, we used the 
values e co = 0.385 kJ/mol and a C o — 3.28 A for the TIP3P, and e C o — 0.395 kJ/mol and 
c"co = 3.26 A for the TIP5P (see [6] for TIP4P parameters), which are similar to those 
derived by Werder et al. [9] to fit the contact angle for a water droplet on a graphene 
surface. A simple analytic form for can be obtained using Steele summation method 
[221 [23] over the graphite periodic structure by writing the interaction of a dispersion 
center with a graphite layer as a Fourier series. The total repulsion-dispersion interaction 
is obtained as a sum of such terms over each graphite layer. We have obtained well 
converged values by including the continuum contribution from the two upper layers and 
the first corrugation from the first layer. 
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In Eq. j2j), V po \ includes the energy associated with the polarization of graphite due to 
the electric field of all the water point charges. This many-body interaction, which turns 
out to be smaller than Vd r , was evaluated using a continuous representation of graphite in 
terms of two contributions, 

V pol = V\\ + V±, (3) 

each one associated, respectively, with the response of graphite to the electric field compo- 
nent parallel and perpendicular to the graphite surface. For the first one, V\\, we assumed 
that graphite behaves as a classical conductor, which allowed us to make use of the image 
charge method to obtain its analytical form. In order to evaluate V±, we associated to 
the graphite surface a constant surface polarizability density a± such that when an electric 
field depending on the surface point and perpendicular to the layer, E±(x, y), is applied, an 
electric dipole density, I(x,y), is induced on that layer, with I(x,y) = a±E±(x,y). Using 
image charge methods one readily shows that if the graphite surface coincides with the 
plane z — 0, the induced image of an electric charge at the point (xi, y^Zi) is an electric 
dipole Pi = —2na±qi at the point [x^y^—z^) and direction parallel to the z axis. This 
result can be generalized additively to the case of several electric point charges to obtain 
an analytical form for V±. The value of the polarizability density a± was estimated from 
£_L, the relative electric permittivity of graphite for applied electric fields perpendicular to 
the (0001) surface, whose value is e± = 5.75; namely, 

where d = 3.35 A is the graphite interlayer distance. We obtained by this procedure 
a± = 0.220 A. 

All other electrostatic contributions to the water-graphite interaction energy having 
vanishing continuous terms (as the water-charge carbon-quadrupole interaction) have been 



neglected, as well as the McLachlan substrate mediated water- water interaction [25]. This 
potential energy surface was argued to be superior to previous empirical models [TO] . 

3 Global Potential Energy Minima 

Likely candidates for the global potential energy minima of graphite-(H 2 0)„ clusters with 
n < 21 were located using the basin-hopping scheme [26], which corresponds to the 'Monte 
Carlo plus energy minimization approach of Li and Scheraga [27J. This method has been 
used successfully for both neutral [26J and charged atomic and molecular clusters [2U 
l28l l29l l30l [31], along with many other applications [32]; of course, this was the method 
used in I. In the size range considered here the global optimization problem is relatively 
straightforward. The global minimum is generally found in fewer than 7x 10 4 basin-hopping 
steps, independently of the random starting geometry. In some cases, starting out from the 
(H 2 0) n global potential minimum, the corresponding global minimum for graphite- (H 2 0) n 
is found even faster. However, the success hit rate of the optimization method decreases 
significantly and the likelihood of our candidates decreases. As a matter of fact, we have 
been able to find for n — 19 and n — 21, in each case, a TIP4P global minimum candidate 
with energy lower than the one found in I; although these new candidates present structures 
very similar to the previously reported ones. 

For graphite- (H 2 0) n clusters, association energies, AE a , are defined for the process 

graphite + nH 2 = graphite-(H 2 0) n ; — (5) 

We also define the water binding energy, AE^, as the difference between the association 
energies of graphite-(H 2 0) n and (H 2 0) n ; that is, 

graphite + (H 2 0)„ = graphite-(H 2 0)„; -AE h . (6) 
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(a) 




(b) 




Figure 1: Two views of the global minimum obtained for graphite- (H 2 0). Water- water 
potential models: TIP3P (a); TIP5P (b). This figure, as well as figures [4] and El was 
prepared using the program XCrysDen |36| . 

The clusters in these expressions are assumed to be in their global minimum. The structures 
and association energies employed here for the global minima of (H 2 0)„ coincide precisely 
with those obtained by Wales and Hodges [TH], Kabrede and Hentschke [33], and James et 
al. PS]. 

In the water monomer case, the structures found for the water-graphite system with 
the TIP3P and TIP5P potentials are given in Fig. Q] (those for TIP4P were presented in 
I). While TIP3P, as well as TIP4P [6], favors a one-legged structure, the TIP5P model 
produces a two-legged global minimum. The equilibrium distance in the global minimum 
between the oxygen and the graphite surface is 3.13 A for the TIP3P potential and 3.12 A 
for the other two models; these distances are very close to the ab initio value (3.04 A) 
[TT] and the corresponding values in water-Ceo (3.19 A) (2TJ and water-benzene (experi- 
mental, 3.33 A) [33] . As happens with the TIP4P model, one-legged and two-legged stable 
structures exist very close in energy for each of the TIPiVP models. Therefore, as we 
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Figure 2: Binding, AE h (full lines), and association, AE^/n (dotted lines), energies 
in kJ/mol for the global minima of water-graphene clusters: TIP3P (triangles), TIP4P 
(squares), TIP5P (stars). 

have discussed in I, the structures found here for the water monomer might change by 
the inclusion in our potential energy surface of the corrugation terms associated with the 
electrostatic interactions, which have been neglected in our PES. For instance, the carbon- 
quadrupole contribution may favor a two-legged structure [34]. However, these effects will 
tend to average out in the adsorption of water clusters. 

The three model potentials provide very similar monomer binding energies, namely 
AE h = 8.81kJ/mol for the TIP4P and AE h = 8.94kJ/mol for the other two models, in 
good agreement with the ab-initio data. The contribution of the polarization energy to 
these binding energies (~ 25%) follows the same trend as the magnitude of the water dipole 
moment for each model and it is responsible for the orientation of the H2O molecule on 
the graphite surface. 

The association (AE a /n) and binding energies (AE b ) for the full graphite-(H 2 0) n clus- 
ters obtained for the three TIPiVP water- water interactions are given in Tableland plotted 
in Fig. El We also show in Fig. [3] the values of the polarization energy V po \ and water- 
graphite dispersion-repulsion energy Vd r , as defined in Section O for the cluster global 
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Table 1: Global minimum association and binding energies in kJ/mol. 
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Figure 3: Dispersion-repulsion, V& r (full line), and Polarization, V^, i (dotted line), contri- 
butions to the binding energies in kJ/mol: TIP3P (triangles), TIP4P (squares), TIP5P 
(stars). 

minima. The three water-water potential models provide very similar values for the as- 
sociation energies. However, the binding energies for TIP4P and TIP5P models differ 
significantly from those obtained with the TIP3P for n > 6. The origin of this difference 
is, as can be seen in Fig. [3J in the Vd r term, which is the dominant contribution for n > 2. 
The term V po \ oscillates with n around an average value of ^ po i = 3.5kJ/mol; the two 
contributions to V po \, V\\ and V±, are similar in magnitude with V\\ somewhat larger than 
V±. The term Vd r fluctuates also around a slowly growing average as the number of water 
molecules close to the graphite surface increases. On average, each of these water molecules 
contributes about 7.3kJ/mol to Vd r - The water-graphite binding energies correspond quite 
closely to the sum of V po \ and Vd r , while the association energies are dominated by the 
water-water interaction. The average value of the association energy per molecule in ho- 
mogeneous TIPiVP (H 2 0) n clusters with 6 < n < 21 is ~ 42kJ/mol [IHICES]. For water 
cluster on graphite the corresponding value turns out to be 44.6kJ/mol, which is compa- 
rable with the experimental value of 43.4±2.9 kJ/mol |7j. Any of these values corresponds 
to the binding energy of a water molecule in a water cluster, and it is much larger than the 
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energy for binding a water molecule onto the graphite surface. This energy balance would 
support an hydrophobic nature of the water-graphite interaction at large scale, as we have 
already discussed in I. 

The structures of the TIP3P lowest minima obtained for graphite-(H20)„ are presented 
in Fig. 01 and those of the TIP5P in Fig. [5] (The corresponding TIP4P global minima 
were presented in I). The three water- water model potentials provide practically identical 
structures for 2 < n < 5. The water substructures in these compounds are actually 
equivalent (see below) to those in the corresponding free global minimum of TIPiVP (H 2 0) n 
[IH [EU [18] , and are in agreement with the ab-initio results. For n — 6, TIP3P and TIP4P 
model have a "book" global minimum, as the one predicted by ab-initio calculations [3lfTT|. 
while the TIP5P leads to an hexagonal ring. This result would favor the first two models 
over the last one. Only the water substructure in the TIP4P differs from that of the 
corresponding free water cluster global minimum ("cage" conformation). 

Here, we will consider two global minima equivalent if their water moieties share the 
same geometrical structure (aside from minor differences in angles and distances) and 
orientation with respect to the graphite surface plane With this convention the equivalences 
found in the global minima of graphite- (H20) n compounds among the three PES have 
been included in the last column of Table CD The number of these equivalences is similar 
to the number of equivalences between the three model potentials in the corresponding 
free water clusters, although these equivalences involve different clusters and potentials. 
For instance, TIP3P and TIP5P provide equivalent free water global minima for n = 8; 
both models present the D 2 d "cube" conformation, while TIP4P global minimum has an 
£4 "cube" conformation. On the other hand, the corresponding water clusters on graphite 
keep the free structure for TIP4P and TIP5P, while TIP3P provides a monolayer water 
conformation; therefore no equivalent structures appear in this case. 

For graphite- water clusters with 6 < n < 17, the TIP3P model seems to favor monolayer 
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water structures (with exceptions for n — 14 and 15), while the other two models favor 
either bilayer (TIP4P and TIP5P) or volume (only TIP5P) structures. The exceptions to 
the monolayer pattern for n = 14, 15 in TIP3P, may be explained by the relative higher 
stability (deduced from second energy differences) of the corresponding free water clusters 
respect to that of their neighbors n = 13, 16. For 17 < n < 21, either volume or bilayer 
conformations are found for the three model potentials. 

When we compare the conformation of the water substructure on the graphite surface 
with that of the corresponding free water cluster, we find also a markedly different behavior 
for the TIP3P model. With the exception of the first six clusters and the case n = 14, those 
two conformations are inequivalent. In other word, the water-graphite interaction is able to 
strongly modify the structure of the free water clusters. This together with the monolayer 
conformation of the adsorbed water clusters would point out to a hydrophilic water-graphite 
interaction for this potential model. However this behavior seem to be due to finite size 
effects since we have not found monolayer global minimum structures for n > 17, neither 
local minima monolayer conformations that are close in energy to the global minima. For 
n > 6 in TIP4P and n > 7 in TIP5P, we do not find either such monolayer structures for 
these two potential models. Thus the hydrophobic nature of the water-graphite interaction 
appears earlier in these models. Furthermore, in some cases for the two latter models, the 
conformation of the water substructure on graphite and that of the corresponding free 
water cluster are equivalent. The exceptions are n — 6, 7, 11, 15, 17, 19, 21 for TIP4P 
and n = 7, n > 13 for TIP5P. In these cases, the water substructure is equivalent to a 
low-lying local minimum of the corresponding TIPiVP (H 2 0) n cluster, rather than to the 
global minimum. The energy penalty for this choice is mainly compensated by a more 
favorable dispersion-repulsion contribution to the interaction energy with graphite, which 
arises from a larger water-graphite contact area. In the structures for the three model 
potentials one finds square and pentagonal water rings; on the other hand, hexagonal rings 
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are less common but they appear more often in TIP5P (n = 11, 12, 14, 15, 19 ) than in 
TIP3P and TIP4P (just for n = 21). 

For TIP4P, the complete two-layer water structures for even n are precisely the struc- 
tures of the global TIP4P free water clusters. Therefore, these structures interact with 
graphite in an optimal way and they keep their structure in the corresponding water- 
graphite clusters. On the other hand, for odd n, the free water global minima do not 
show optimal surfaces for its interaction with graphite, thus explaining why these clusters 
change their structure to minimize that interaction energy. The chosen new structures are 
sensibly determined by those of either the n — 1 or n + 1 clusters. The TIP5P (and also 
TIP3P) potential model do not produce this alternating behavior in the structure of the 
free water global minima and, therefore, we find a different behavior in the water-graphite 
global minima for n > 8. 

Second energy differences account for the relative cluster stability; their values for 
association and binding energies, per water molecule, are plotted in Fig. El TIP4P and 
TIP5P, show practically the same behavior in the whole n range. The n = 4 cluster is 
particularly stable in all cases. For n > 10, we observe that the three model potentials 
present an oscillation of period An = 2, namely, clusters with even n are more stable than 
their odd n neighbors. This is an interesting feature because it does not occur so neatly for 
the free water clusters, and it is not obviously related, except in the TIP4P case [6], with 
the cluster structures. Second differences for the binding energies do not show common 
patterns for the three potential models because these provide very different global minima 
structures respect to the water-graphite interaction. 
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Figure 6: Second energy differences per water molecule (in kJ/mol) for the association 
energies (a) and binding energies (b) of water-graphite clusters. Water-water potential 
models: TIP3P (triangles, full line), TIP4P (squares, dotted line), TIP5P (stars, dashed 
line) 
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4 Conclusions 



Using basin-hopping global optimization and a potential energy surface built up from three 
different water- water interaction models (TIP3P, TIP4P and TIP5P) we have characterized 
the geometrical structures and energetics of the likely candidates for the global potential 
energy minima of graphite- (H 2 0) n clusters up to n = 21. The structures of these minima 
for 1 < n < 5 coincide for the three potential models with those provided by other available 
calculations. The global minimum for the compound with n = 6 agrees with the ab-initio 
structure for TIP3P and TIP4P, but not for TIP5P. For n > 6, no ab-initio data are 
available and, except for the equivalences presented in Table I, the three model potentials 
provide different global minimum structures, as occurs for the free water clusters. For 
n > 2, association energies are dominated by the water- water interaction while the main 
contribution to the binding energies comes from the dispersion energy; furthermore the 
polarization term V po \ can be safely neglected for the larger clusters (n > 3); this justifies 
the use of water-graphite potentials that include only dispersion-repulsion terms [9j. For 
small n, the water grows on the graphite surface forming a monolayer. However, as n 
increases the hydrophobic nature of the water-graphite interaction dominates and breaks 
this tendency. The threshold for this transition is at n = 7 for TIP4P and TIP5P and 
n = 18 (with the exceptions n = 14, 15) for TIP3P. Therefore this latter potential seem to 
favor planar conformations up to larger n. 

The hydrophobic character of the water-graphite interaction at the nanoscopic level 
makes in some cases the water substructure in the lowest energy clusters to be equivalent to 
a low-lying minimum of the appropriate (H20) n free cluster. In many cases the structure is 
simply a slightly relaxed version of the global minimum for (H20) n , and therefore equivalent 
to it. For TIP3P this occurs only for the first six n values and for n = 14. TIP5P shows 
equivalences for the same first six n values and for n — 8, 9, 10, 11, 12. TIP4P shows the 
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larger number of equivalences for n < 5 and n — 8, 9, 10, 12, 13, 14, 16, 18, 20. 

The lowest energy structures obtained in the present work will be made available for 
download from the Cambridge Cluster Database |3oJ. 
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